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We study the Bose-Einstein condensate trapped in a three-dimensional spherically symmetrical 
potential. Exact solutions to the stationary Gross-Pitaevskii equation are obtained for properly 
modulated radial nonlinearity. The solutions contain vortices with different winding numbers and 
exhibit the shell-soliton feature in the radial distributions. 
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The experimental realization of matter- wave solitons in weakly interacting atomic Bose-Einstein condensates (BECs) 
has attracted great attention in the past decade. Dark[l], bright [2], and gap solitons[3], as well as vortices [4, l5| have 
been observed and widely studied. By applying the Feschbach resonance techniques, scientists can conveniently control 
the nonlinear it ies. Many nonlinear phenomena have been predicted by the manipulation of the scattering length either 
in time[7,Q or iri^ space j9|,|l3|. It has been used to generate bright solitons[10, 11] and induce collapse[12]. In addition, 
bright solitons [14] and periodic wave solutions [15] were obtained in spinor BECs. However, thus far, only effectively 
one-dimensional matter-wave solitons have been created in experiments which used cigar-shaped traps to confine 
the condensate [l6|. The realization of higher-dimensional matter- wave solitons in BECs is still a challengeable topic 
because those solutions are usually unstable for the nonlinear Schrodinger (NLS) equation with constant or uniform 
couplings due to the weak and strong collapses [6]. Modulation of atomic scattering length by the Feshbach resonance 
is expected to dynamically stabilize higher dimensional bright solitons p^. 

On the other hand, a multiquantum vortex is typically dynamically unstable in harmonically trapped BEC |18M2QL 
It will split into single quantum vorticesjiue to the splitting instability [21, 22.]. In the presence of a plug potential [23|, 
an anharmonic trapping potential 0, [25|], or a immiscible two-species BECQ, the multiquantum vortex can be 
effectively stabilized. For example, when the confining potential is steeper than harmonic potential, multiquantum 
vortices are energetically favorable. 

In this paper we study 3D solutions to the Gross-Pitaevskii equation (GPE) by means of the similar transformation. 
Exact solitary vortices in the BEC trapped in an external potential are constructed by properly modulating the 
spatially modulated nonlinearity. The number of vortex-soliton (VS) modes is specified by the radial nodes in the 
wavefunction and the winding number of the vortices. 

The scaled stationary GPE for the macroscopic wave function reads 

M^(r) = -V2^(r) + 5(r)|^(r)| V(r) + F(r)^(r), (1) 

where V(v) is the spherically symmetric well which is chosen as the harmonic oscillator potential, ^(r) is the nonlin- 
earity coefficient which is spatially variable. By rewriting the wavefunction in the spherical coordinates as 

i,{r,9,^) = <^{r)Yi^{e,^), (2) 

where Yim{0^ (p) = P^'^' (cos6>)e*^'^ (/ = 0, 1, 2, • • • , and m = 0, ±1, ±2, • • • , ±1) is the spherical harmonic function, the 
radial part ^{r) obeys the following nonlinear equation 

^$ = _ + g(r)$3 + [fcil + y(r)]$. (3) 

The converge condition for r ^ requires that ^(r) ~ for / ^ 0, while ^\r) = for / = 0. In the meantime, 
$(r ^ oo) due to the localization of the wavefunction. 

According to the similar transformation, we define ^{r) = p{r)U[R{r)] so as to obtain two independent 
equations pJi], 
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and 



goU^ = EU, 



(5) 



where R{r) = x '^{x)dx which can be viewed as a rescaled radius. For this purpose, the spatial modulation 
of nonlinear coupling should be chosen as g{r) = gor~^p~^{r). Here E and go are constants. The nonlinearity- 
modulation is determined by p(r). Equation (j5j) is exactly solvable in terms of Jacobian elliptic functions which gives 
rise to solitary form of wave. For Eq.(|4]), physically meaningful solutions impose restrictions on p(r). As r ^ 0, p{r) 
behaves as with a > 2/3. It should diverge for r ^ oo so that the nonlinearity g{r) is bounded and the integration 
in R{r) converges. 

We first consider the repulsive nonlinearity (^o > 0) by including the harmonic trap V{r) = ar'^ to confine the 
system. The existence of solutions to Eq.([5j) requires E > 0. The VS solution is 



U{R) = ^/2{E-k^)/gosn{kR,p), 



(6) 



where k is the period of the Jacobian elliptical function sn. The modulus p = y^^^T^^^^ imposes restric- 
tion on k as \/E/2 < k < \fE. To meet boundary conditions <^>(0) = <l>(oo) = 0, it requires that k = 
nK{^E/k'^ — l/R{oo)) (n = 2,4,6, •••). Here K{p) is the complete elliptic integral of the first kind. It follows 
from \/E /2 < nK{^E/k'^ — l/R{oo)) < \fE that the number of VS n < nmax = 2R{oo)^/E / n ^ which implies that 
there is only a finite number of radial VS modes (or none, if nmax < 2). 

On the other hand, Eq.(|l]) is a nonlinear equation which can be solved numerically. To perform the calculation, 
one requires p ~ r"'^' at r ^ 0. Thus the nonlinear term in Eq.(6) may be neglected near r = 0. In this case, p 
behaviors like the Neumann function Ys{^/Jir) at r ^ for /i > (for /i < it can be checked that the VS solutions 
do not exist). On the other end of r ^ oo, p ^ oo due to the presence of the harmonic trap. In addition, Er~^p~^ 
with ^ > in Eq.(j4]) guarantees the sign definiteness of p{r). Numerically, one takes the initial distribution of p{r) 
by setting a small value of ro to the Neumann function Ys{-sJJlr) and its derivative. The function R{r) and g{r) can 
be simultaneously obtained by adjusting the value of s according to the angular momentum /. 




FIG. 1: (a) The rescaled radius R{r) for / = 1 (blue line) and / = 2 (red line), respectively, (b) Modulation profiles of the 
repulsive nonlinearity g{r) for / = 1 (blue line) and / = 2 (red line), respectively. The parameters are ^ = 1, /x = 2, a = 0.01, 
and go = 0.01. 



Figure 1 display the numerical results of radial distributions of the rescaled radius R{r) and the modulated non- 
linearity g{r) for £^ = 1, /i = 2, a = 0.01, and go = 0.01. We have chosen s = 1.38 for / = 1 and s = 2.56 for / = 2, 
respectively. It can be seen from that R{r) is angular momentum-dependent and saturates as r ^ oo. In Fig. 1(b), 
the nonlinear modulation exhibits oscillations versus the radius. It vanishes at both r = and large distances. In 
addition, g{r) is dependent on the angular momentum Figure 2 shows the analytical result for U{R). Obviously, 
it satisfies the restrictions U{R{0)) = U{R{oo)) = for both / = 1 (Fig2.(a)) and / = 2 (Fig.2(b)). The number of 
nodes increases as the radial quantum number n increases. The radial wavefunction ^n^(^) = Pi{'^)Uni{r) are shown 
in Fig. 3. One notes that the radial quantum number n should be even numbers. 

The radial wavefunction ^(r) is plotted in Fig. 3 for / = 1 (Fig. 3 (a)) and / = 2 (Fig. 3(b)), respectively. We see it is 
localized in the trap. We observe that the number of nodes is exactly equal to that of U{r) as shown in Fig. 2. This 
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FIG. 2: U{R) versus r for (a) / = 1, and (b) / = 2. Open circles (blue), squares (red), solid circles (black), and triangles (green) 
denote for n — 2, 4, 6, 8, respectively. The parameters are the same as in Fig.l. The number of nodes equals to n/2 — 1. 




FIG. 3: The radial function $(r) for (a) / = 1 and (b) / = 2 in the presence of the harmonic trap. Open circles (blue), squares 
(red), solid circles (black), and triangles (green) denote for n = 2, 4, 6, 8, respectively. The parameters are the same as in Fig.l. 



means that the function p{r) does not contribute nodes to the radial wavefunction. The whole 3D wavefunction is the 
combination of the radial and the angular part of the wavefunction, tpnimi'^^^^^) = ^ni(^)^Zm(^7 It is specified 
by three quantum numbers (n/2,/,m). The 3D density distributions are shown in Fig. 4, where the upper raw and 
lower raw are for the radial quantum numbers n = 2 and n = 6, respectively. Fig. 4 (a) and (b) are for the angular 
momentum / = 1 (m = 0, 1, respectively). Fig. 4(c)- (e) are for / = 2 (m = 0, 1, 3, respectively). 

For attractive nonlinearity (^o < 0), the results are similar to those of the repulsive nonlinearity as £^ 7^ 0. 
Here we pay our attention to the case of £^ = 0. Since the nonlinear term disappears, Eq.(|3]) is exactly solvable 
for the harmonic trap V{r) = ar^. p{r) can be represented in terms of the Whittaker' M and W Functions as 
p{r) = r-^/2[ciM(/i/4yS, (|/ + l/2|)/2, y^r^) + C2W{p/4^/a, {\l + l/2|)/2, y/ar^)]. In the absence of the external 
potential, p{r) degenerates into p{r) = c^ji{\/ — pr) + C4nl{^/—pr) with /i < 0, where ji and ni are the spherical Bessel 
and spherical Neuman functions of imaginary argument. In the meantime, Eq.(j5j) for £^ = yields 

U{R) = {nr]/^^)cn{nr]R - K{l/V2), 1/V^), (7) 

where the radial quantum number n = 2, 4, 6, • • • and r] = K {\/2 / 2) / R{oo) . 

The results in the presence of the harmonic potential are displayed in Fig. 5. The relevant parameters are chosen 
as/i = — 1, ^o = — 1, and a = 0.01. The number of nodes in the radial wavefunctions are determined by the number 
of nodes in U (r). 

In summary, we have studied exact solitary solutions to the spherical symmetric GPE with radially modulated 
nonlinearity. The wavefunction can be expressed by the combination of a radial function and a spherical harmonic 
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FIG. 4: The 3D density distributions of the states ([2]) for n = 2 (upper raw) and n = 6 (lower raw), (a) / = l,m = 0; (b) 
/ = 1, m = 1; (c) / = 2, m = 0; (d) / = 2, m = 1; (e) / = 2, m = 2. The parameters are the same as in Fig.l. 
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FIG. 5: Results for the attractive nonlinearity with ^ = and / = 1. (a) The rescaled radius R{r). (b) The modulated 
nonlinearity g{r). (c) U{r) versus r. (d) The radial distribution of the wavefunction. Open circles (blue), squares (red), solid 
circles (black), and triangles (green) denote for n = 2,4,6,8, respectively. 



function. The number of VS are limited and the quantum number can only be even numbers. 

This work is supported by the funds from the Ministry of Science and Technology of China under Grant No. 
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